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Abstract 

We investigate dynamic properties of a plasma whose constituents are confining 
gluons resulting from the Gribov-Zwanziger quantisation. In a static formula¬ 
tion, this system reproduces qualitatively the pure-glue equation of state and 
thereby encodes crucial features of the phase transition. The dynamic descrip¬ 
tion proposed in this work allows us to study close-to-equilibrium transport 
phenomena with the inclusion of confinement effects. In particular, we deter¬ 
mine the non-equilibrium behaviour of the interaction measure (trace anomaly) 
and find the form of the bulk viscosity coefficient. The latter can be used in 
phenomenological applications to heavy-ion collisions. 


1. Introduction 

Recent experimental and theoretical studies of ultrarelativistic heavy-ion 
collisions suggest that the quark-gluon plasma formed in these processes be¬ 
haves like a strongly-interacting and dissipative fluid [1]. These observations 
have triggered a broad interest in the relativistic hydrodynamics of viscous sys¬ 
tems and their respective transport coefficients [2]. The latter, by force of the 
Kubo formulas, are sensitive to the long-distance dynamics of the underlying 
microscopic theory [3]. The expected running-coupling strength in the phe¬ 
nomenologically relevant regime is however on the order of unity, g ~ 0{1) [4], 
it is therefore not a surprise that many of the experimental features are unsat¬ 
isfactorily described by conventional perturbation theory, see [5] for the most 
up-to-date results on thermodynamic quantities. Moreover, the dynamic na¬ 
ture of the plasma evolution severely limits the applicability of lattice methods 
[6]. Hence the understanding of deconfinement remains a pressing challenge for 
heavy-ion physics. 
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A way to tackle these obstacles is to improve the description of soft non- 
Abelian gauge degrees of freedom. In this work we introduce for the first time 
a non-equilibrium and dynamic description of a plasma consisting of confining 
gluons obtained from the Gribov-Zwanziger quantization of Yang-Mills (YM) 
theory [7]. In this scenario, fixing the infrared (IR) residual gauge transforma¬ 
tion in the Coulomb gauge generates a new scale that leads to an IR improved 
dispersion relation for gluons [7], 



( 1 ) 


where k is the gluon three-momentum and E the gluon energy. The Gribov 
parameter 7 g, which governs the onset of non-Abelian effects in the IR, is found 
self-consistently from a gap equation [7, 8]. This leads to a suppression of soft 
degrees of freedom which is a necessary condition for color conhnement [7, 8, 
9]. Thus low-momentum gluons are confining while their ultraviolet properties 
remain unaltered [8]. 

The Gribov framework has been extensively studied in the vacuum, see 
Ref. [10] for reviews, where it is in good agreement with lattice and functional 
methods, especially in the Coulomb gauge [11, 12]. The generic modification of 
the IR properties of YM, embodied in Eq. (1), is however operative at any tem¬ 
perature, and the framework has also been generalized to finite temperature in 
order to study the deconfined phase, see e.g. [13, 14, 15, 16]. In the phenomeno¬ 
logically relevant temperature regime the Gribov parameter 7 g is approximately 
a constant [15] and, recently, its possible relation with the phase transition has 
also been investigated [17]. The resulting equation of state provides a good 
description of lattice data down to the vicinity of the critical temperature 
[13, 15]. In particular, the interaction measure, see Eq. (5), exhibits a char¬ 
acteristic peak near the phase transition [13]. It has also been realized that 
this mechanism induces long-range correlations through novel massless collec¬ 
tive modes exclusively existing for YM theory [16, 18]. Given these advances, 
it is intriguing to explore a dynamical application of the Gribov framework, 
especially bearing in mind the pressing challenges from heavy-ion experiments, 
and in the following we provide a first effort in this direction. 

We will denote a hot system of gluons, governed locally by Eq. (1), as the 
Gribov-Zwanziger (GZ) plasma. Based on the discussion above, we argue that 
essential confinement effects are intrinsically incorporated in our setup once 
the Gribov dispersion relation is adopted. Since our present study focuses on 
the phenomenologically relevant regime, we will set 7 g = const, motivated by 
Refs. [13, 15]. 

Eirstly, we consider local thermodynamic properties of the GZ plasma and, 
secondly, determine its dynamic evolution in a boost-invariant and transversally 
homogeneous (0 -I- 1)D system [19]. This allows us to find exact solutions for 
in- and out-of-equilibrium evolutions using the kinetic theory approach in the 
relaxation time approximation [20, 21, 22]. 

One of our main results is to determine the connection between the Gribov 
parameter 7 g and the bulk viscosity C,. While the ratio of bulk viscosity to 
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entropy density C/s is found to be strongly suppressed for gluons at asymp¬ 
totically high temperatures [23], it is expected to be enhanced in the vicinity 
of the YM phase transition [24], see also [25, 26] for results in quasi-particle 
models. This is in line with lattice results [27] and a significant bulk viscosity 
is also obtained from calculations in strongly-coupled theories [28]. Recently 
the influence of bulk viscosity on observables in heavy-ion collisions [29, 30], as 
well as on the stability of the hydrodynamic evolution [31], has attracted a lot 
of attention. Our results highlight for the first time the direct connection of ( 
to the confining properties of YM theory in the vicinity of the phase transition 
and its significance for phenomenological applications to heavy-ion collisions. 


2. Covariant kinetic equation setup 

In order to generalize the results [13] to a dynamical formulation, one needs 
to implement Lorentz covariance while respecting the dispersion relation Eq. (1) 
in the local rest frame. The dispersion relation in Eq. (1) may hence be rewritten 
in covariant form, 

E{k ■ u) = ^{k ■ uY + ')‘^/{k ■ u)2 , (2) 

where u is the four-velocity of the fluid element.^ The use of Eq. (2) allows us 
to restrict our considerations to energies that are invariably defined in the local 
rest frame by Eq. (1) where the Gribov parameter is evaluated. Further details 
about our covariant setup can be found in [32]. Thus, generalizing the results of 
[13], the energy density and pressure of the fluid described by the distribution 
function / = f{x, k) are given in covariant forms by 


e = 


dK E{k ■ u) f, 


( 3 ) 


P = 


1 

3 


dK 


{k ■ uY 
E{k ■ u) 



/, 


( 4 ) 


respectively, where the integration measure is dK = god^k{k ■ u)/[(27r)^fc°] and 
go = 2{N^ — 1) is the degeneracy factor for gluons with Nc colors {go = 16 
for SU(3)). These expressions can be used to find the interaction measure 
(sometimes referred to as the trace anomaly). 


I = e-3P = 2 


dK 


7g 


(fc • u)2 E{k ■ u) 


/• 


( 5 ) 


^We assume that the in-medium value of 70 is determined in the fluid element’s local 
rest frame, where = (1,0, 0,0). We introduce = |fc|, which is the magnitude of the 
three-vector k = (kx,ky,k^^), and kj_ = (k‘^ -h such that the resulting four-vector 

k^ = (/c®, k) has standard Lorentz transformation properties with = 0. The four-vector k 
can be interpreted as a four-momentum of a perturbative, non-interacting gluons. Equation 
(2) seems to be a straightforward generalisation of the non-covariant expression (1), however, 
other version of (1) might be also possible. 
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This quantity is closely related to the deviation from conformality of the system 
under consideration. We note that I vanishes identically in the 7 g 0 limit, 
i. e., for a conformal gas of massless quasi-particles. 

In the case of local equilibrium, the function / has the Bose-Einstein form 
/gz = [exp(E(fc-u)/T(a:)) —1]“^, where the temperature T can depend on space 
and time. The energy density and pressure obtained from Eqs. (3) and (4) with 
/ = /gz will be denoted below as Egz and Pgz, respectively. The temperature 
dependence may in this case be eliminated to construct the equation of state of 
the perfect GZ plasma, Egz = £gz(Tgz)- 

In order to study dynamic phenomena away from equilibrium, we employ 
kinetic theory in the relaxation time approximation (RTA) [21]. The (0 -I- 1)D 
system is described by the Bjorken flow vector = x^/t = (t/r, 0, 0, z/r) 
where r = is the proper time. Eurthermore, we introduce the boost- 

invariant variables [33] v = kPt — k\\z and w = feyt — kPz, where k-u = 
vjr = \J Consequently, the integration measure becomes AK = 
go dr(;d^A:_L/[(27r)^r]. We note that the thermodynamic functions now depend 
only on the proper time, while the distribution function may depend on r, w 
and kx- 

Calculating the proper time derivative of the energy density given by Eq. (3), 
which in these coordinates reads s{t) = f dK E{t, w, k±) /(r, w, k±), one gets 


de £ -I- P|| 
^ 

dr T 


= J dKE{T,w,kx) 


df{T,w,k±) 

dr 


( 6 ) 


where we identify the component of the pressure acting in the longitudinal 
direction as 


P|| = / dK 


w 


t'^E{t, w, k±) 


1 - 


7^ 


{up- E k\Y 


/• 
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With w‘^/{t‘^E) replaced here by k‘j_/{2E) one may find the transverse pressure. 


P± 



kj 

2E{t, w, k±) 


1 - 


7^ 


{up -I- 


/, 


( 8 ) 


and check that P — {2P± -|- P||)/3.^ 

The terms on the left-hand side of Eq. (6) cancel due to the energy-momentum 
conservation, = 0. This implies that the term on the right-hand side in 


® For a (0 -I- 1)D system, the energy-momentum tensor is 

= (£ + P_l)uI^u'' - + (P|| - Px)z^P, 

where = (t, 0, 0, z)/r and = {z, 0, 0, t)/T. The energy density and the two pressures are 
given by the formulas e = 3P = 2Px + ^|| = —P|| = z^z^^T^^^ where 

. In our case, the energy density and the pressures are expressed by the 
integrals over the distribution function /, see Eqs. (3), (7), and (8), respectively. For more 
details see [20, 21, 22] . 
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Eq. (6) should vanish as well and suggests the use of the standard RTA kinetic 
equation of the form 

df{T,w,k±) ^ fGz{r,w,k±) - f{T,w,k±) 
dT ~ rrei(r) ’ 

altogether with the Landau matching condition, Eqz = £, given explicitly by 
f dKE(T,w,k±) {fez — f) =0. We emphasise that Eq. (9) is applicable for 
the close-to-equilibrium situations but this is sufficient for derivation of our 
central result, namely, the formula for the bulk viscosity.^ 

The formal solution of Eq. (9) is [20, 21, 22] 

r dr' 

f{T,W,k±) = fo{w,k±)D{T,To) + J ■^^^-^D{T,T')fGz{T',W,k±), (10) 

where the distribution function at the initial proper time tq is given by 
/(tq, w, k±) = /o(w, k±) and the damping function is given by 


D{t 2 ,ti) = exp 





( 11 ) 


In order to construct the solution of Eq. (10) we have to know the dependence 
of T and Tj-ei on the proper time r. Although the relaxation time Trei can, in 
general, depend on the temperature, in this work we fix it to a constant in order 
to single out the genuinely novel features arising in the GZ plasma. 

In the Trei —0 limit, the form of Eq. (9) guarantees that the actual distribu¬ 
tion function tends rapidly to the equilibrium one and Py = P± = Pqz- Thus, in 
order to obtain the leading-order temperature profile, we solve the well-known 
Bjorken hydrodynamic equation 


d£oz(P(r)) £az{T{T)) + Poz{T{t)) 
dr r 


which follows directly from the energy-momentum conservation [19]. In this 
case, the equilibrium entropy density is Sqz = (egz + Pgz)/T, where deoz = 
Tdsez and dP^z = SezdT. The temperature dependence of the perfect GZ 
plasma, which we denote Tgzj follows directly from these relations, and is found 
by solving {Tgz{t)), where cf = i9Pgz/9£gz is the speed of 

sound of the plasma. For conformal systems, where e = 3P, we reproduce 
the well-known scaling solution T(r) = To{tq/tY= with Cg = 1/3. In the GZ 
plasma, Cg tend to the ideal value only in the high-P limit while deviating from 
it in the studied temperature range. 

Away from equilibrium the temperature of the system is determined from 
the Landau matching condition by demanding that / yields the same energy 


^Free-streaming in Eq. (9) is recovered by taking the vanishing coupling limit, correspond¬ 
ing to simultaneously taking the relaxation time to infinity and the Gribov scale to zero, see 
e.g. [15]. 
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Figure 1: Evolution of the temperature of the system. Depicted are the results for perfect 
fluid (black, thin line), equilibrium GZ plasma (GZP) (red, solid) and two non-equilibrium 
evolutions characterized by the relaxation times Trei =1 fm/c (dashed, green curve) and Trei = 
2 fm/c (dotted, blue curve) 


density as the GZ equilibrium function Taking the appropriate moment of 
the solution of the kinetic equation (10) we find the integral equation 


Egz (T(r)) = n(r,To)i?e(^,^ 




Here we have introduced the auxiliary functions 

yhe{y,h) 


Heia,b) = 


9olt 


dy 


27r2 Jq exp[a^/y^~TT/y^] - 1 


(14) 


and h^{y,b) = b /y^d(/)sin(0) + 1 with I3'^{b,(j)) = 6^ cos^ </) + sin^ (j)- 

The temperature dependence can then implicitly be read off from the left-hand 
side of Eq. (13) by taking advantage of the equilibrium relation Egz (T) = 
The evolution of the pressure is found in a completely analo¬ 
gous way, i.e., by inserting the solution given in Eq. (10) into the right-hand 
side of Eq. (4) with the temperature dependence found in the previous step. 


3. Thermodynamics-like quantities 

We proceed with numerical calculations based on Eq. (10), where the energy 
density and pressure are read off from Eqs. (3) and (4), respectively. Eor the 
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Figure 2: Temperature dependence of the interaction measure I scaled by T^. Upper panel: 
The result for the GZ plasma (GZP) in local equilibrium (solid, red curve) together with the 
lattice data [34] (black dots). Lower panel: The two non-equilibrium results for I, character¬ 
ized by the relaxation times r^ei =1 fm/c (dashed, green curve) and r^ei = 2 fm/c (dotted, 
blue curve), compared to the equilibrium calculation (solid, red curve). 


• /QCD (Borsanyi et al, JHEP 07, 056 (2012)) 


7 g = 700 MeV 

• \ 

• 

• 

- GZP (eq.) 


(a) : 

• ' 11 

- 

.w*, . 




initial condition we choose the GZ plasma in equilibrium at a given temperature. 
As stated in the discussion above, we fix the Gribov parameter to the value 
7 g = 700 MeV in order to obtain the best description of lattice data at T > Tc 
in agreement with previous studies [13]. Furthermore, we choose the initial 
time To = 0.5 fm/c, the initial temperature T{to) = 600 MeV, and study the 
non-equilibrium evolutions for two fixed relaxation times, Trei = 1 fm/c and 
Trel = 2 fm/c. 

The proper time dependence of the temperature is shown in Fig. 1. The GZ 
plasma (GZP) in equilibrium (solid, red curve) cools at a slower rate than the 
ideal gas (black, thin curve) owing to the reduction of the speed of sound in the 
vicinity of the critical temperature. The rate of cooling for local equilibrium is 
only mildly restrained when going away from equilibrium, see the dashed, green 
and dotted, blue curves in Fig. 1 which correspond to the two equilibration times 
described above, respectively (we keep the same color coding in all subsequent 
figures). 

The upper panel of Fig. 2 shows our results for the interaction measure, 
Eq. (5), compared to the lattice data [34], where we have chosen = 260 MeV. 
The solid red line depicts the equilibrium result, which qualitatively describes 
the lattice data: it yields half of the peak of the anomaly in the region of the 
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phase transition and gives a good description of the anomaly in the temperature 
range just above the phase transition, T ~ 1.5 — ATc- In this way we reproduce 
the result of Ref. [13] . Further improvements can only be achieved by taking into 
account quantum effects [15]. The dashed and dotted lines in the lower panel of 
Fig. 2 depict the results of the non-equilibrium calculations for the two chosen 
values of Trei and approach the solid, equilibrium line at lower temperatures, 
reflecting that for late times r ^ Tj-ei the system approaches equilibrium. The 
results presented in the lower panel of Fig. 2 indicate that the trace anomaly 
becomes larger if the system evolves out of equilibrium. In addition, the width 
of the peak becomes wider. Similarly to the equilibrium results, we expect that 
our kinetic results are the most legitimate in the same temperature region. 


4. Bulk viscosity from the kinetic equation 

We note that the interaction measure vanishes exactly at all temperatures 
for a conformal theory. In our case, it is the presence of an IR non-Abelian scale 
in the dispersion relation Eq. (1), which allows us to qualitatively reproduce the 
YM theory around the phase transition. Although the energy density is kept 
equal to the equilibrium one due to the Landau matching condition, the pres¬ 
sure can deviate. The difference between the actual and equilibrium pressures 
describes the bulk viscous pressure 11 = P — P^z ■ Owing to the definition of the 
trace anomaly in Eq. (5) and the Landau matching, one finds 

which simply corresponds to (/qz — The knowledge of 11 allows us to 

determine the effective bulk viscosity ^eff of the system through the identihcation 


n = -Ceff = -Ceff/r , 


(16) 


where the latter equality holds for a (0 -I- 1)D system. 

For small deviations from equilibrium, the effective viscosity (^eff approaches 
the standard bulk viscosity coefficient C. To determine C we seek the solution of 
the kinetic equation (9) in the form / « /gz + Sf. In the linear approximation, 
we find Sf = —Treid/cz/dr and substitute Sf into the right-hand side of Eq. (15). 
Using the equilibrium relation for the speed of sound, we find 


c 


go7G Pel 

37r2 T 



-1g 

+ 7l 


/gz(1 + /gz), 


(17) 


where /gz = [exp(-\/g^ + ~ l] which is one of the main results 

of our paper. The bulk viscosity is proportional to the relaxation time by 
construction. It vanishes when 7 g —t 0, which is the case for a conformal, 
massless gas. Our setup is hence qualitatively different from the outset due to 
the presence of the Gribov parameter. 









Figure 3: Left panel: The effective bulk viscosity evaluated from the kinetic theory, Eqs. (15) 
and (16) (solid curves), and the bulk viscosity coefficient obtained from the linearized ex¬ 
pression in Eq. (17) (dashed and dotted curves)). Right panel: Temperature dependence of 
the bulk viscosity coefficient scaled by the equilibrium entropy density for two values of the 
relaxation time. The color coding is the same as in Figs. 1 and 2. 


In the left panel of Fig. 3 we show the proper time dependence of the ef¬ 
fective bulk viscosity Ceff obtained directly from the numerical solution of the 
kinetic equation using Eqs. (15) and (16), and of the bulk viscosity coefficient C 
obtained from Eq. (17). As the system approaches equilibrium we find a good 
agreement between and C: which supports the validity of Eq. (17). More¬ 
over, it is worth noting the close relationship between C and the deviation from 
conformality, understood here as cf — 1/3, see Eq. (17). On a qualitative level, 
this is also observed in quasi-particle models [26], however the magnitude of 
C/^Gz is significantly smaller in this case. Further improvements are necessary 
for the GZ gas in order to facilitate more qualitative studies. In a parallel work 
[32], we have derived a scaling of the ratio of bulk to shear viscosities to the 
speed of sound, C/rj = «:gz(1/3 — c^) at high temperature that is similar to the 
results from calculations in holographic, strongly-coupled theories [28]. 

In the right panel of Fig. 3, on the other hand, we show the ratio C/^gz as 
a function of temperature. Similarly to Fig. 1, we observe a significant increase 
of the bulk viscosity in the vicinity of the phase transition. We note that large 
bulk viscosity may be related to prolonged deviations from equilibrium [29]. 
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5. Summary 

In this work, we have established a novel framework for dealing with a non¬ 
equilibrium plasma of confining gluons in close to equilibrium situations. Our 
method is based on the improved dispersion relation for gluons, Eq. (1), which 
accounts for residual long-range correlations in the system in the deconfined 
phase. If local equilibrium is assumed, our system expands and cools as ex¬ 
pected from the Bjorken model, see the red curve in Fig. 1. In non-equilibrium 
situations, our numerical results show an increase of the trace anomaly at a 
given effective temperature, see Fig. 2. In the latter case, the system evolu¬ 
tion is affected by dissipative phenomena, in particular, by the presence of bulk 
viscosity. The full and linearized solutions, derived in Eq. (17), for the bulk 
viscosity agree at large times and indicate a rise of the ratio C,/s close to the 
phase transition as expected in a strongly coupled plasma. 

The features above provide an improved understanding of the hot plasma 
produced in heavy-ion collisions, in line with expectations from lattice data, 
and implies that bulk viscosity should be implemented in their hydrodynamic 
modeling. Finally, we note that the framework proposed herein can be naturally 
extended to include other transport phenomena such as, for example, the shear 
viscosity, and be systematically improved in order to account for a more involved 
and realistic form of the relaxation time and a temperature dependence of the 
Gribov parameter [14, 15]. We leave these developments, which would allow for 
more qualitative assessments of the dynamical properties of the GZ plasma, for 
future works. 
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